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We study statistical scale invariance and dynamic scaling in a simple solid-on-solid 2 + 1— dimen- 
sional limited mobility discrete model of nonequilibrium surface growth, which we believe should 
describe the low temperature kinetic roughening properties of molecular beam epitaxy. The model 
exhibits long-lived transient anomalous and multiaffine dynamic scaling properties similar to that 
found in the corresponding 1 + 1— dimensional problem. Using large-scale simulations we obtain 
the relevant scaling exponents, and compare with continuum theories. 



A key issue in kinetic surface roughening (l],|^] is mak- 
ing connection between theoretical growth universality 
classes, as defined by continuum growth equations for ex- 
ample, and experimentally observed rough growth in real 
surfaces which generally depends on many details such as 
growth conditions (eg. temperature, surface orientation, 
growth rate) and atomistic rules controlling attachment, 
detachment, evaporation, and most importantly, diffu- 
sion of the deposited adatoms at the growth front. The 
concept driving much of the kinetic surface roughening 
research is that only a few universality classes de- 
scribe the asymptotic growth properties in many seem- 
ingly different nonequilibrium surface growth problems 
as most of the details are irrelevant from a renormaliza- 
tion group viewpoint and do not affect the asymptotic be- 
havior. Much recent work has gone into building simple 
atomistic discrete nonequilibrium growth models which 
catch the essential aspects of a complicated growth prob- 
lem and include only the relevant dynamical processes 
determining the asymptotic growth behavior. One such 
nonequilibrium growth model was introduced by one of 
us in rcf. || in the context of one dimensional molecular 
beam epitaxy. This growth model has since been exten- 
sively studied and it seems to describe Q well the 
low temperature growth properties of realistic stochastic 
Monte Carlo simulation results of molecular beam epi- 
taxy. Although the growth model introduced in ref. @] 
has been fairly extensively studied in the literature M, 
almost all of the existing work is in 1 + 1 dimensions 
where a one dimensional substrate roughens as it grows. 
We present in this paper results of a systematic study 
of the growth model of ref. |^| in the physically relevant 
2 + 1 dimensions. 

In our growth model ||] , atoms are randomly deposited 
on an I x I (we have studied system sizes upto L = 10 
with a maximum of 10 7 deposited layers, which amount 
to the deposition of upto 10 13 atoms) flat substrate un- 
der solid-on-solid deposition and growth conditions. If a 
randomly deposited atom has at least one lateral nearest- 
neighbor bond (i.e. if its initial coordination number is 
2 or more), then it is incorporated at the deposition site 
and stays there forever. Otherwise the atom could move 
to a nearest-neighbor lateral site (with no restriction on 
the number of vertical sites it moves in the growth di- 



rection) for incorporation provided it can increase (but 
not necessarily maximize) its coordination number at the 
final site. If no such nearest-neighbor lateral site is avail- 
able (which would increase the atom's coordination num- 
ber) the atom is again incorporated at the site of de- 
position. If more than one final site could increase its 
coordination number, the atom randomly moves to any 
one of these final sites with equal probability and is in- 
corporated there permanently. The motivation for this 
manifestly nonequilibrium growth model is that during 
low temperature molecular beam epitaxy it is unlikely 
for atoms to break lateral bonds, and only deposited 
adatoms without any lateral bonds can move with ap- 
preciable mobility. This is a typical example of a limited 
mobility growth model where atoms at kink and trapping 
sites (i.e. those with at least one lateral bond) simply do 
not move. A typical example of a saturated growth mor- 
phology resulting from these rules is shown in Fig.l. 




FIG. 1. The saturated growth morphology on a 100 x 100 
substrate at 3 x 10 6 monolayers. Darker(lighter) shades rep- 
resent lower(higher) points on the surface. 

We first study the global scaling behavior of the surface 
by considering the dynamical surface width W , which is 
the root-mean-square height fluctuation of the growing 
surface, defined as 

W(L,t) = {(h- (h))Y /2 - (1) 

If there is statistical scale invariance in the problem, then 
W scales with the lateral system size L and time t as 
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W(L,t)~L<f(L/at)), 
where the scaling function f(y) is 



m = 



const for i/<Cl 
for y ^> 1. 
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Here £ is the roughness exponent, and £(t) is the lateral 
correlation length obeying the dynamic scaling relation 
~ i 1 / 2 where z is the dynamical exponent. Time 
in the simulation is measured in the average number of 
deposited layers. Combining these scaling relations the 
interface width W can be written as W(L, t <C L z ) ~ t 13 
and W(L,t ^> L z ) ~ L^, where /3 = C/ z 1S the growth 
exponent. Our calculated results for W(L, t) are shown 
in Fig. 2 for a substrate size L = 500 and 100 (inset). The 
calculated growth exponent (for L = 500) clearly shows a 
crossover from an initial value of about 0.25 to an asymp- 
totic value of about 0.2. We find the same crossover be- 
havior in other large systems (L — 500-1000) we have 
studied — in smaller systems {L = 100), however, the 
crossover to the asymptotic exponent (/3 « 0.2) is not 
seen and j3 remains around 0.25 as can be seen in the 
inset of Fig. 2. In a second inset, the values of the satu- 
ration width W(L,t — -> oo) are plotted as a function of 
the system size L = 20 to 100. (The small L values in 
the saturation plot reflect the high value of z ~ 3 in the 
problem : L = 100 requires more than 10 6 layers to sat- 
urate.) The slope yields the global roughness exponent 
C = 0.56 ±0.1. 
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FIG. 2. The interface width W as a function of deposition 
time t in the 500 x 500 system (left inset: 100 x 100 system). 
Right inset: The saturation width W sa t vs the substrate size 
L. 

The local scaling behavior is studied by calculating 
the height-height correlation function G(r, t) defined as 
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G(r,t) = (\h(x + T,t) - h(x,t)\ 2 ) 
scaling form of the correlation function is 

G(r,t)~r</(r/tft)) 



The conventional 
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where r — |r| and the asymptotic behavior of the scal- 
ing function /(y) is the same as that of the function / 



in Eq.(3). It is, however, known that the correspond- 
ing one dimensional growth model exhibits anomalous || 
and multiafhne scaling behavior in contrast to the sim- 
ple conventional dynamic scaling defined in Eq.(3). We 
find the same to be true in the 2+1 dimensional model 
also with the only difference being that the anomalous 
behavior is manifestly a long-lived (over three decades in 
time) transient in 2 + 1 dimensions whereas in 1 + 1 di- 
mensions it lasts for at least eight decades in time 
(and perhaps much longer; the true asymptotic limit may 
not have yet been reached in 1 + 1 dimensions Qj) . In or- 
der to study multifractality, we follow Krug || and define 
a generalized correlation function by considering higher 
moments of the equal-time height difference correlator 



G q (r,t) = (|/i(x + r,t)-/i(x,i)| 9 > 1/9 - 



(5) 



These functions G q (note that for q = 2 one gets back 
the G of Eq.(4)) have g-dependent roughness exponents 
if the growing surface is multiaffine whereas they all scale 
with the same exponent in a selfaffine statistically scale 
invariant surface. 

In Fig. 3 we show the multifractal behavior of G q (r,t) 
in our model. We find exactly the same qualitative and 
quantitative multifractal behavior (with the same criti- 
cal exponents within numerical accuracy) for the different 
system sizes (L = 100 — 1000) we study. We summarize 
below (see ref. for details) the asymptotic behav- 

ior of the correlation function G q in the anomalous and 
multiaffine scaling situation depicted in Fig. 3. 

Multiaffine scaling (as opposed to statistically scale in- 
variant selfaffine scaling) signifies ^-dependent dynamic 
scaling properties of G q : 



, . r*« for r < £(t) 
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In addition, the local and the global dynamic scaling prop- 
erties are necessarily different in this anomalous scaling 
situation. Thus, / in Eq.(4) scales differently from / in 
Eqs.(2) and (3): 



f(y) 



for y < 1 
for y > 1. 
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Using the asymptotic behavior of £(t) we obtain the 
asymptotic behavior for the gth moment of the height- 
height correlation function as 



G q (r,t) 




for r < t 1/z < L 
for r <C L < t 1/z 
for r > t 1/z . 



(8) 



Note that an additional scaling exponent a q is needed to 
fully characterize the anomalous multiscaling situation. 
In the conventional scaling situation a q = and Q q = Q. 

In Fig. 3(a) we show our results for G q for L = 500 at 
10 3 monolayers. It is obvious that the lines for different q 
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FIG. 3. (a) Height-height correlation function G q (r), 
q — 1 — 4 from bottom to top. (L = 500 at 10 3 monolayers). 
The solid lines indicate power law fits with slope C,' q as shown 
in the table, (b) Nearest neighbor height difference function 
a g (t), q = 1 — 4 from bottom to top (L = 1000). Inset: The 
local derivative of logcr^ with respect to \ogt. O, A, □ and x 
correspond to q = 1, 2, 3 and 4 respectively. 



have different local roughness exponents, C,' q = C,q~ a q, im- 
plying multifractality. The local roughness exponents (' q 
vary from 0.184 for q = 4 to 0.395 for q = 1. Comparing 
C 2 = 0.3 with C = 0.56 (cf. Fig. 2) we conclude that the lo- 
cal and global roughness in our problem have different be- 
havior (anomalous scaling Q), which is one direct conse- 
quence of the spatial multiscaling. In order to gain more 
insight into the multiaffine scaling behavior we show in 
Fig. 3(b) the time evolution of the nearest-neighbor height 
difference and its higher moments a q (t) = G q (r = l,i) 
for L = 1000. Substituting r = 1 in Eq.(8), we get (fol- 
lowing ref. |g|) 



t ot q /z for t l/z 

L a " for t 1 ! 2 



< L 
> L. 



(9) 



In the conventional scaling situation a q = for all q, 
and therefore a q (t) quickly saturates to a constant value 
almost immediately. But in the anomalous scaling case 
a q (t) grows with time before saturating (at very 

long time) in the steady state when £(t) ~ L. The ap- 
proximate crossover time for this saturation is t c ~ L z 
and can be very large [||||||]. In Fig. 3(b) the nearest 



neighbor height difference a q increases with time for all 
q in the small time region. For 1 < t < 20, the slope 
(a q /z) ranges from 0.22 (q=l) to 0.31 (q=4). It is abso- 
lutely clear, however, from Fig. 3(b) that o~ q (t) is showing 
a crossover behavior where the effective exponent (a q /z) 
of Eq.(9) is approaching zero long before the steady-state 
saturation regime where £(t) ~ L. It is obvious from the 
inset that the slopes of a q (t) curves, instead of staying at 
constant values of a q /z, are going down to zero around 
t ~ 10 5 , which is long before the expected saturation time 
t c > 10 9 (L=1000). (The corresponding a q (t) results for 
a smaller L — 100 system size show the vanishing of the 
effective exponents a q /z around t « 10 3 , implying that 
multiscaling lasts for about three decades in time in a 
100 x 100 system.) The bending of the cr q (t) curve as a 
function of i, which is glaringly apparent in our Fig.3(b), 
can also be seen (in a much less pronounced fashion) in 
the corresponding 1 + 1 dimensional results Our 
results presented in Fig. 3(b) make it obvious that the 
multiscaling observed in our 2+1 dimensional simula- 
tions is a long-lived transient and is not the asymptotic 
behavior of the limited mobility epitaxial growth model 
introduced in ref. 0. Our results strongly indicate that 
the same must be true in 1 + 1 dimensions also, which has 
already been suggested in ref. M, except that the mul- 
tiscaling transient in 1 + 1 dimensions lasts for at least 
10 8 time steps and perhaps much longer. We have, in 
fact, studied a 1 + 1 dimensional version of our prob- 
lem for L = 10 upto 10 8 monolayers of growth, and we 
find that a q / z decreases from an initial value of a q /z = 
0.195( 9 = l),0.222(g = 2),0.262(g = 3),0.294( 9 = 4) to 
a q j z — .107, .143, .168, .179 definitively establishing that 
the multiscaling behavior ^-^| of the growth model in- 
troduced in ref. || is an extremely long-lived transient, 
even in 1 + 1 dimensions. 

Having established that the anomalous multiscaling 
behavior of our limited mobility growth model, while be- 
ing extremely interesting, is most definitely a transient 
(and not an asymptotic) behavior, we need to ask the 
obvious question : What is the universality class of the 
discrete growth model introduced in ref. (^]? The asymp- 
totic critical exponents (3 ~ 0.2 and ( « 0.6 in Fig. 2 are 
consistent with the nonlinear MBE growth equation in- 
troduced in refs. || and |JlC L We believe that the dis- 
crete growth model of ref7~| 1 , which we study here in 
2 + 1 dimensions, does indeed asymptotically belong to 
the universality class of the nonlinear growth equation 

ifiol 



^^ 4 V 4 ft + A 2 V 2 (V/i) 2 
at 



(10) 



which in 2 + 1 dimensions has [|10| the critical exponents 
(3 = 1/5 and £ = 2/3, which, within numerical errors, 
is consistent with the critical exponents we obtain in 
this paper pl[ . The corresponding linear equation (with 
A2 = in Eq.(10)) has (3 = 1/4 in 2 + 1 dimensions, ex- 
plaining the crossover from (3 w 0.25 to 0.20 seen in our 
Fig. 2. (It should be emphasized that in 1 + 1 dimensions 
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even the largest simulations of this model ^-0] obtain 
[3 w 0.375 which is consistent with the linear growth 
equation (A2 = 0), and the anticipated crossover to the 
nonlinear model has not yet been definitely established.) 
Note that Eq.(10) breaks the up-down symmetry in the 
problem and is therefore a true noncquilibrium model of 

growth. The measured skewness 5* 



- {(h-(h))*)*/* 



and 



the effective fourth cumalant Q 



3 in our 



simulations are approximately —0.5 and 1.2 respectively 
in the saturated regime. 

The anomalous scaling and multifractality in our re- 
sults most likely arise from the existence of an infinite 
series of higher-order marginally irrelevant terms of the 



on the right hand side of Eq.(10), 



form^A 2 „V 2 (V/i) 2 ™ 

n=2 

as has recently been speculated in ref. Q and discussed 
extensively in ref. [Q . Further discussion JjJ of this issue 
is beyond the scope of our work, and we only note that 
this infinite series of higher order gradient terms should 
have a much stronger influence in 1 + 1 dimensions where 
they are all marginally relevant, explaining why the ob- 
served multifractality is so much more pronounced in 1+1 
dimensions ||. 

Finally, we note that the standard Laplacian term, 
V2S7 2 h, of the Edwards- Wilkinson equation Jl2| is absent 
for the growth model introduced in ref. |5| in contrast to 
other popular solid-on-solid growth models Jll||l4| stud- 
ied in the literature in the context of epitaxial growth. 
The absence of the V 2 /i term (and its fourth order coun- 
terpart JTot ] , the V(V/i) 3 term) in the growth model of 
ref. H is, in fact, an exact result arising from a hid- 
den symmetry [0 in the problem which makes the in- 
clination dependent current on a tilted substrate [l5| in 
this growth model vanish exactly (our calculated current 
on a tilted substrate in our 2 + 1 dimensional simula- 
tion is ±10~ 4 , which is zero within our numerical ac- 
curacy). We have also studied the closely related Wolf- 
Villian model JlJ], where the deposited adatoms (inde- 
pendent of their initial lateral coordination) seek out the 
nearest-neighbor sites of maximum coordination, which 
is now definitely known [jl5 16 1 to asymptotically belong 
to the Edwards- Wilkinson universality class, where our 
measured tilt dependent surface current has a small neg- 
ative (downhill) value (increasing in magnitude with in- 
clination) of approximately —0.1 at a slope of 2. We see 
essentially no multiscaling behavior in our 2 + 1 dimen- 
sional Wolf- Villain model simulations : a q / z for L = 500 
is essentially q independent and decreases from a very 
small initial value of about .06 to zero at t « 10 3 . Clearly, 
the Wolf- Villain model behaves very differently from the 
model of ref. (3) in 2 + 1 dimensions even though their 
effective behavior is known to be 'almost' identical in 
1 + 1 dimensions! The fact that these two models, ref. 
H and pH |, have different universality classes was actu- 
ally pointed out some time ago [[l7| although substan- 
tial confusion still seems to exist about their relation- 



ship. Results presented here strongly support the idea 
that the limited mobility noncquilibrium epitaxial growth 
model introduced in ref. || is asymptotically described 
by the nonlinear MBE growth equation (Eq.10) although 
there are interesting non-asymptotic corrections associ- 
ated with the anomalous multiscaling phenomena. The 
significance of this finding lies in the fact that most exper- 
imental investigations J_8| of the kinetic surface roughen- 
ing phenomenon in epitaxial growth obtain large growth 
((3 k 0.2 — 0.3) and roughness (a » 0.5 — 1) exponents, 
indicating that the linear and the nonlinear MBE growth 
equation may be playing important roles in real growth. 
The model of ref. Q seems to be the only known solid-on- 
solid epitaxial growth model which is asymptotically not 
described by the linear second-order Edwards- Wilkinson 
equation, but by the fourth-order nonlinear MBE growth 
equation. 
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